import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
modelling_home = '/Volumes/mo2016/home/Documents/modelling'
modelling_local = root + '/Documents/modelling'
import matplotlib as mpl
mpl.use('tkagg')
if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
modelling_ephemeral = root + '/ephemeral/Documents/modelling'
modelling_home = root + '/home/Documents/modelling'
modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
print(root)
sys.path.append(modulepath)
from plotting_numerical import *
from tqdm import tqdm
import pickle
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import matplotlib.animation as animation
import matplotlib as mpl
#############################
#Opening list with parID's
lhs=True
if lhs==True:
folder = 'fullcircuit_5716gaussian/1M_turingI'
variant=0
details = '1M_turingI'
else:
var=float(sys.argv[1])
# var=0.23
folder = 'fullcircuit_5716gaussian/var%s'%var
variant='5716gaussian'
details = 'var%r'%var
circuit_n=2
shape='ca'
mechanism = 'fullcircuit'
dimension='2D'
L=10; x_gridpoints =15; J = L*x_gridpoints
T =120; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/ca/%s'%(folder)
if lhs==True:
if folder == 'fullcircuit_5716gaussian/1M_turingI':
print('yes')
general_filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,'1MturingI', shape,mechanism,L,J,T,N)
else:
general_filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
general_filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)
# metric=str(sys.argv[1])
metric='ps_min'
parID_dict = pickle.load( open(modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/%s_dict_%s.pkl"%(metric,general_filename), "rb" ) )
print(modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/%s_dict_%s.pkl"%(metric,general_filename))
len(parID_dict)
parID_list = []
entropy_list = []
parID_dict = dict(sorted(parID_dict.items(), key=lambda item: item[1])) #important to sort out dictionary by HKS values
for key in parID_dict:
parID_list.append(key)
entropy_list.append(parID_dict[key])
# print(parID_list, entropy_list)
parID_list = [int(i) for i in parID_list] #turn string list into integer list
print(parID_list)
# k=20
num=len(parID_list)
n_col = int(np.sqrt(num))
# for i in range(0,k):
# n_col = 10
# row = np.where(fit==i)[0]
# print(row) # row in Z for elements of cluster i
# num = row.shape[0] # number of elements for each cluster
n_row = np.floor(num/n_col)+1 # number of rows in the figure of the cluster
# print("cluster "+str(i))
# print(str(num)+" elements")
fig = plt.figure(figsize=(n_col/10+2, n_row/10+2))
dx = float(L)/float(J-1)
grid = np.array([j*dx for j in range(J)])
for count,parID in tqdm(enumerate(parID_list[:10]),disable=False):
ax=plt.subplot(n_row,n_col, count+1)
#rgb_timeseries=timeseries_unstacked_list[row[n]] # Read the numpy matrix with images in the rows
# par_ID = parID_list[row[n]]
if lhs==True:
filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
else:
filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
final_concentration = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
rgb = plot_redgreenblue_contrast(final_concentration,L,mechanism,shape,filename,modelling_ephemeral,parID=parID,scale_factor=x_gridpoints,save_figure='LargeImage')
# rgb_timeseries=timeseries_unstacked # Read the numpy matrix with images in the rows
# ax.set_title(parID,size=0.1)
ax.set(yticks=[],xticks=[],yticklabels=[],xticklabels=[])
ax.imshow(rgb.astype('uint8'), origin= 'lower')
ax.set_ylabel('%s-%s'%(parID,entropy_list[count]),size= 1,c='y', labelpad=0.35)
# plt.title('1M numerical search 0-%r'%num)
# plt.savefig(modelling_home + '/3954/numerical_confocal/results/entropy/LargeImages/%s_%s.png'%(metric,general_filename), dpi=2000)
plt.show()
# plt.savefig('h.png',dpi=2000)
print('gh')
# plt.clf()
# plt.close(fig)
Running cells with 'Python 3.7.3 64-bit' requires ipykernel package.
Run the following command to install 'ipykernel' into the Python environment.
Command: '/usr/bin/python3 -m pip install ipykernel -U --user --force-reinstall'
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
modelling_home = '/Volumes/mo2016/home/Documents/modelling'
modelling_local = root + '/Documents/modelling'
import matplotlib as mpl
mpl.use('tkagg')
if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
modelling_ephemeral = root + '/ephemeral/Documents/modelling'
modelling_home = root + '/home/Documents/modelling'
modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)
from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################
# %matplotlib inline
def fourier2d(final_concentration):
# print('Normalised PowerSpectrum of original')
pixels=final_concentration[5]
pixels[pixels == 0] = 0.0000000001
pixels/=np.sum(pixels) #line added
#fourier transform
f = np.fft.ifftshift(pixels)
f = np.fft.fft2(f)
f = np.fft.fftshift(f)
if plot==True:
plt.subplot(142)
plt.imshow(np.real(f))
plt.colorbar()
plt.subplot(143)
plt.imshow(np.real(f))
plt.colorbar()
# max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
PowerSpectrum = np.real(f)**2 + np.imag(f)**2
# PowerSpectrum[max_coordinates]=0
if plot==True:
plt.subplot(144)
plt.imshow(PowerSpectrum)
plt.colorbar()
if np.sum(PowerSpectrum)!=0:
PowerSpectrum/=np.sum(PowerSpectrum)
# print('sum',np.sum(PowerSpectrum))
# print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))
# plt.subplot(132)
# plt.imshow(PowerSpectrum)#,norm=LogNorm())
# plt.colorbar()
# binwidth = 0.00001
PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
# print(np.sum(PowerSpectrum_vector))
# plt.subplot(133)
# hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
# plt.yscale('log')
# plt.show()
return PowerSpectrum_vector
def entropy(histogram):
x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
return x
lhs=True
if lhs==True:
folder = 'fullcircuit_5716gaussian/1M_turingI'
variant=0
else:
var=float(sys.argv[1])
# var=0.23
folder = 'fullcircuit_5716gaussian/var%s'%var
variant='5716gaussian'
circuit_n=2
shape='square'
mechanism = 'fullcircuit'
L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )
parID_ps = {}
plot=False
# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
# for parID in tqdm(parID_list[:20], disable=False):
for parID in tqdm([805,686,472,252,688], disable=True):
# print('parID',parID)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
else:
filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
if not np.isnan(final_concentration0).any():
final_concentration = np.round(final_concentration0,4)
if plot==True:
plt.figure(figsize=[14,2])
plt.subplot(141)
plt.imshow(final_concentration[5])
plt.colorbar()
hist_PowerSpectrum = fourier2d(final_concentration)
ps = entropy(hist_PowerSpectrum)
parID_ps[parID] = ps
# print(final_concentration0)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)
pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )
# if root=='/rds/general/user/mo2016':
# sendmail('general_filename')
Running cells with 'Python 3.7.3 64-bit' requires ipykernel package.
Run the following command to install 'ipykernel' into the Python environment.
Command: '/usr/bin/python3 -m pip install ipykernel -U --user --force-reinstall'
Python 3.8.5 (default, Sep 4 2020, 07:30:14)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.19.0 -- An enhanced Interactive Python. Type '?' for help.
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
modelling_home = '/Volumes/mo2016/home/Documents/modelling'
modelling_local = root + '/Documents/modelling'
import matplotlib as mpl
mpl.use('tkagg')
if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
modelling_ephemeral = root + '/ephemeral/Documents/modelling'
modelling_home = root + '/home/Documents/modelling'
modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)
from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################
# %matplotlib inline
def fourier2d(final_concentration):
# print('Normalised PowerSpectrum of original')
pixels=final_concentration[5]
pixels[pixels == 0] = 0.0000000001
pixels/=np.sum(pixels) #line added
#fourier transform
f = np.fft.ifftshift(pixels)
f = np.fft.fft2(f)
f = np.fft.fftshift(f)
if plot==True:
plt.subplot(142)
plt.imshow(np.real(f))
plt.colorbar()
plt.subplot(143)
plt.imshow(np.real(f))
plt.colorbar()
# max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
PowerSpectrum = np.real(f)**2 + np.imag(f)**2
# PowerSpectrum[max_coordinates]=0
if plot==True:
plt.subplot(144)
plt.imshow(PowerSpectrum)
plt.colorbar()
if np.sum(PowerSpectrum)!=0:
PowerSpectrum/=np.sum(PowerSpectrum)
# print('sum',np.sum(PowerSpectrum))
# print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))
# plt.subplot(132)
# plt.imshow(PowerSpectrum)#,norm=LogNorm())
# plt.colorbar()
# binwidth = 0.00001
PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
# print(np.sum(PowerSpectrum_vector))
# plt.subplot(133)
# hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
# plt.yscale('log')
# plt.show()
return PowerSpectrum_vector
def entropy(histogram):
x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
return x
lhs=True
if lhs==True:
folder = 'fullcircuit_5716gaussian/1M_turingI'
variant=0
else:
var=float(sys.argv[1])
# var=0.23
folder = 'fullcircuit_5716gaussian/var%s'%var
variant='5716gaussian'
circuit_n=2
shape='square'
mechanism = 'fullcircuit'
L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )
parID_ps = {}
plot=False
# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
# for parID in tqdm(parID_list[:20], disable=False):
for parID in tqdm([805,686,472,252,688], disable=True):
# print('parID',parID)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
else:
filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
if not np.isnan(final_concentration0).any():
final_concentration = np.round(final_concentration0,4)
if plot==True:
plt.figure(figsize=[14,2])
plt.subplot(141)
plt.imshow(final_concentration[5])
plt.colorbar()
hist_PowerSpectrum = fourier2d(final_concentration)
ps = entropy(hist_PowerSpectrum)
parID_ps[parID] = ps
# print(final_concentration0)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)
pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )
# if root=='/rds/general/user/mo2016':
# sendmail('general_filename')
/rds/general/user/mo2016
--------------------------------------------------------------------------- FileNotFoundError Traceback (most recent call last) ~/Documents/modelling/3954/numerical_confocal/code/entropy/plot_all_colony.py in <module> 110 else: 111 filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N) --> 112 final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) ) 113 if not np.isnan(final_concentration0).any(): 114 final_concentration = np.round(final_concentration0,4) FileNotFoundError: [Errno 2] No such file or directory: '/rds/general/user/mo2016/home/Documents/modelling/3954/numerical_confocal/results/simulation/square/fullcircuit_5716gaussian/1M_turingI/2Dfinal_circuit2_variant0_square_fullcircuitID805_L5_J50_T2000_N20000.pkl'
#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
modelling_home = '/Volumes/mo2016/home/Documents/modelling'
modelling_local = root + '/Documents/modelling'
import matplotlib as mpl
mpl.use('tkagg')
if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
modelling_ephemeral = root + '/ephemeral/Documents/modelling'
modelling_home = root + '/home/Documents/modelling'
modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)
from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################
%matplotlib inline
# %matplotlib inline
def fourier2d(final_concentration):
# print('Normalised PowerSpectrum of original')
pixels=final_concentration[5]
pixels[pixels == 0] = 0.0000000001
pixels/=np.sum(pixels) #line added
#fourier transform
f = np.fft.ifftshift(pixels)
f = np.fft.fft2(f)
f = np.fft.fftshift(f)
if plot==True:
plt.subplot(142)
plt.imshow(np.real(f))
plt.colorbar()
plt.subplot(143)
plt.imshow(np.real(f))
plt.colorbar()
# max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
PowerSpectrum = np.real(f)**2 + np.imag(f)**2
# PowerSpectrum[max_coordinates]=0
if plot==True:
plt.subplot(144)
plt.imshow(PowerSpectrum)
plt.colorbar()
if np.sum(PowerSpectrum)!=0:
PowerSpectrum/=np.sum(PowerSpectrum)
# print('sum',np.sum(PowerSpectrum))
# print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))
# plt.subplot(132)
# plt.imshow(PowerSpectrum)#,norm=LogNorm())
# plt.colorbar()
# binwidth = 0.00001
PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
# print(np.sum(PowerSpectrum_vector))
# plt.subplot(133)
# hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
# plt.yscale('log')
# plt.show()
return PowerSpectrum_vector
def entropy(histogram):
x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
return x
lhs=True
if lhs==True:
folder = 'fullcircuit_5716gaussian/1M_turingI'
variant=0
else:
var=float(sys.argv[1])
# var=0.23
folder = 'fullcircuit_5716gaussian/var%s'%var
variant='5716gaussian'
circuit_n=2
shape='square'
mechanism = 'fullcircuit'
L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )
parID_ps = {}
plot=False
# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
# print('parID',parID)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
else:
filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
if not np.isnan(final_concentration0).any():
final_concentration = np.round(final_concentration0,4)
if plot==True:
plt.figure(figsize=[14,2])
plt.subplot(141)
plt.imshow(final_concentration[5])
plt.colorbar()
hist_PowerSpectrum = fourier2d(final_concentration)
ps = entropy(hist_PowerSpectrum)
parID_ps[parID] = ps
# print(final_concentration0)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)
pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )
# if root=='/rds/general/user/mo2016':
# sendmail('general_filename')
15%|█▌ | 3/20 [00:00<00:00, 23.29it/s]
/rds/general/user/mo2016
100%|██████████| 20/20 [00:00<00:00, 22.90it/s]
#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
modelling_home = '/Volumes/mo2016/home/Documents/modelling'
modelling_local = root + '/Documents/modelling'
import matplotlib as mpl
mpl.use('tkagg')
if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
modelling_ephemeral = root + '/ephemeral/Documents/modelling'
modelling_home = root + '/home/Documents/modelling'
modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)
from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################
%matplotlib inline
plot=True
# %matplotlib inline
def fourier2d(final_concentration):
# print('Normalised PowerSpectrum of original')
pixels=final_concentration[5]
pixels[pixels == 0] = 0.0000000001
pixels/=np.sum(pixels) #line added
#fourier transform
f = np.fft.ifftshift(pixels)
f = np.fft.fft2(f)
f = np.fft.fftshift(f)
if plot==True:
plt.subplot(142)
plt.imshow(np.real(f))
plt.colorbar()
plt.subplot(143)
plt.imshow(np.real(f))
plt.colorbar()
# max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
PowerSpectrum = np.real(f)**2 + np.imag(f)**2
# PowerSpectrum[max_coordinates]=0
if plot==True:
plt.subplot(144)
plt.imshow(PowerSpectrum)
plt.colorbar()
if np.sum(PowerSpectrum)!=0:
PowerSpectrum/=np.sum(PowerSpectrum)
# print('sum',np.sum(PowerSpectrum))
# print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))
# plt.subplot(132)
# plt.imshow(PowerSpectrum)#,norm=LogNorm())
# plt.colorbar()
# binwidth = 0.00001
PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
# print(np.sum(PowerSpectrum_vector))
# plt.subplot(133)
# hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
# plt.yscale('log')
# plt.show()
return PowerSpectrum_vector
def entropy(histogram):
x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
return x
lhs=True
if lhs==True:
folder = 'fullcircuit_5716gaussian/1M_turingI'
variant=0
else:
var=float(sys.argv[1])
# var=0.23
folder = 'fullcircuit_5716gaussian/var%s'%var
variant='5716gaussian'
circuit_n=2
shape='square'
mechanism = 'fullcircuit'
L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )
parID_ps = {}
# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
# print('parID',parID)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
else:
filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
if not np.isnan(final_concentration0).any():
final_concentration = np.round(final_concentration0,4)
if plot==True:
plt.figure(figsize=[14,2])
plt.subplot(141)
plt.imshow(final_concentration[5])
plt.colorbar()
hist_PowerSpectrum = fourier2d(final_concentration)
ps = entropy(hist_PowerSpectrum)
parID_ps[parID] = ps
# print(final_concentration0)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)
pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )
# if root=='/rds/general/user/mo2016':
# sendmail('general_filename')
5%|▌ | 1/20 [00:00<00:02, 7.11it/s]
/rds/general/user/mo2016
100%|██████████| 20/20 [00:02<00:00, 7.06it/s]
#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
modelling_home = '/Volumes/mo2016/home/Documents/modelling'
modelling_local = root + '/Documents/modelling'
import matplotlib as mpl
mpl.use('tkagg')
if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
modelling_ephemeral = root + '/ephemeral/Documents/modelling'
modelling_home = root + '/home/Documents/modelling'
modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)
from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################
%matplotlib inline
plot=True
# %matplotlib inline
def fourier2d(final_concentration):
# print('Normalised PowerSpectrum of original')
pixels=final_concentration[5]
pixels[pixels == 0] = 0.0000000001
pixels/=np.sum(pixels) #line added
#fourier transform
f = np.fft.ifftshift(pixels)
f = np.fft.fft2(f)
f = np.fft.fftshift(f)
if plot==True:
plt.subplot(142)
plt.imshow(np.real(f))
plt.colorbar()
plt.subplot(143)
plt.imshow(np.real(f))
plt.colorbar()
max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
PowerSpectrum = np.real(f)**2 + np.imag(f)**2
PowerSpectrum[max_coordinates]=0
if plot==True:
plt.subplot(144)
plt.imshow(PowerSpectrum)
plt.colorbar()
if np.sum(PowerSpectrum)!=0:
PowerSpectrum/=np.sum(PowerSpectrum)
# print('sum',np.sum(PowerSpectrum))
# print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))
# plt.subplot(132)
# plt.imshow(PowerSpectrum)#,norm=LogNorm())
# plt.colorbar()
# binwidth = 0.00001
PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
# print(np.sum(PowerSpectrum_vector))
# plt.subplot(133)
# hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
# plt.yscale('log')
# plt.show()
return PowerSpectrum_vector
def entropy(histogram):
x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
return x
lhs=True
if lhs==True:
folder = 'fullcircuit_5716gaussian/1M_turingI'
variant=0
else:
var=float(sys.argv[1])
# var=0.23
folder = 'fullcircuit_5716gaussian/var%s'%var
variant='5716gaussian'
circuit_n=2
shape='square'
mechanism = 'fullcircuit'
L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )
parID_ps = {}
# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
# print('parID',parID)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
else:
filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
if not np.isnan(final_concentration0).any():
final_concentration = np.round(final_concentration0,4)
if plot==True:
plt.figure(figsize=[14,2])
plt.subplot(141)
plt.imshow(final_concentration[5])
plt.colorbar()
hist_PowerSpectrum = fourier2d(final_concentration)
ps = entropy(hist_PowerSpectrum)
parID_ps[parID] = ps
# print(final_concentration0)
# if lhs==True:
# filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
# else:
# filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)
pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )
# if root=='/rds/general/user/mo2016':
# sendmail('general_filename')
5%|▌ | 1/20 [00:00<00:01, 9.63it/s]
/rds/general/user/mo2016
100%|██████████| 20/20 [00:02<00:00, 8.28it/s]
#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
modelling_home = '/Volumes/mo2016/home/Documents/modelling'
modelling_local = root + '/Documents/modelling'
import matplotlib as mpl
mpl.use('tkagg')
if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
modelling_ephemeral = root + '/ephemeral/Documents/modelling'
modelling_home = root + '/home/Documents/modelling'
modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)
from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################
%matplotlib inline
plot=True
save=True
if plot==True:
save==False
# %matplotlib inline
def fourier2d(final_concentration):
# print('Normalised PowerSpectrum of original')
pixels=final_concentration[5]
pixels[pixels == 0] = 0.0000000001
pixels/=np.sum(pixels) #line added
#fourier transform
f = np.fft.ifftshift(pixels)
f = np.fft.fft2(f)
f = np.fft.fftshift(f)
if plot==True:
plt.subplot(142)
plt.imshow(np.real(f))
plt.colorbar()
plt.subplot(143)
plt.imshow(np.real(f))
plt.colorbar()
max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
PowerSpectrum = np.real(f)**2 + np.imag(f)**2
PowerSpectrum[max_coordinates]=0
if plot==True:
plt.subplot(144)
plt.imshow(PowerSpectrum)
plt.colorbar()
if np.sum(PowerSpectrum)!=0:
PowerSpectrum/=np.sum(PowerSpectrum)
# print('sum',np.sum(PowerSpectrum))
# print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))
# plt.subplot(132)
# plt.imshow(PowerSpectrum)#,norm=LogNorm())
# plt.colorbar()
# binwidth = 0.00001
PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
# print(np.sum(PowerSpectrum_vector))
# plt.subplot(133)
# hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
# plt.yscale('log')
# plt.show()
return PowerSpectrum_vector
def entropy(histogram):
x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
return x
lhs=True
if lhs==True:
folder = 'fullcircuit_5716gaussian/1M_turingI'
variant=0
else:
var=float(sys.argv[1])
# var=0.23
folder = 'fullcircuit_5716gaussian/var%s'%var
variant='5716gaussian'
circuit_n=2
shape='square'
mechanism = 'fullcircuit'
L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )
parID_ps = {}
# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
# print('parID',parID)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
else:
filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
if not np.isnan(final_concentration0).any():
final_concentration = np.round(final_concentration0,4)
if plot==True:
plt.figure(figsize=[24,4])
plt.subplot(141)
plt.imshow(final_concentration[5])
plt.colorbar()
hist_PowerSpectrum = fourier2d(final_concentration)
ps = entropy(hist_PowerSpectrum)
parID_ps[parID] = ps
# print(final_concentration0)
if save==True:
if lhs==True:
filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)
pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )
# if root=='/rds/general/user/mo2016':
# sendmail('general_filename')
5%|▌ | 1/20 [00:00<00:02, 7.31it/s]
/rds/general/user/mo2016
100%|██████████| 20/20 [00:02<00:00, 7.77it/s]
#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
modelling_home = '/Volumes/mo2016/home/Documents/modelling'
modelling_local = root + '/Documents/modelling'
import matplotlib as mpl
mpl.use('tkagg')
if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
modelling_ephemeral = root + '/ephemeral/Documents/modelling'
modelling_home = root + '/home/Documents/modelling'
modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)
from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################
%matplotlib inline
plot=True
save=True
if plot==True:
save==False
# %matplotlib inline
def fourier2d(final_concentration):
# print('Normalised PowerSpectrum of original')
pixels=final_concentration[5]
pixels[pixels == 0] = 0.0000000001
pixels/=np.sum(pixels) #line added
#fourier transform
f = np.fft.ifftshift(pixels)
f = np.fft.fft2(f)
f = np.fft.fftshift(f)
if plot==True:
plt.subplot(142)
plt.imshow(np.real(f))
plt.colorbar()
plt.subplot(143)
plt.imshow(np.imag(f))
plt.colorbar()
max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
PowerSpectrum = np.real(f)**2 + np.imag(f)**2
PowerSpectrum[max_coordinates]=0
if plot==True:
plt.subplot(144)
plt.imshow(PowerSpectrum)
plt.colorbar()
if np.sum(PowerSpectrum)!=0:
PowerSpectrum/=np.sum(PowerSpectrum)
# print('sum',np.sum(PowerSpectrum))
# print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))
# plt.subplot(132)
# plt.imshow(PowerSpectrum)#,norm=LogNorm())
# plt.colorbar()
# binwidth = 0.00001
PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
# print(np.sum(PowerSpectrum_vector))
# plt.subplot(133)
# hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
# plt.yscale('log')
# plt.show()
return PowerSpectrum_vector
def entropy(histogram):
x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
return x
lhs=True
if lhs==True:
folder = 'fullcircuit_5716gaussian/1M_turingI'
variant=0
else:
var=float(sys.argv[1])
# var=0.23
folder = 'fullcircuit_5716gaussian/var%s'%var
variant='5716gaussian'
circuit_n=2
shape='square'
mechanism = 'fullcircuit'
L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )
parID_ps = {}
# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
# print('parID',parID)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
else:
filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
if not np.isnan(final_concentration0).any():
final_concentration = np.round(final_concentration0,4)
if plot==True:
plt.figure(figsize=[24,4])
plt.subplot(141)
plt.imshow(final_concentration[5])
plt.colorbar()
hist_PowerSpectrum = fourier2d(final_concentration)
ps = entropy(hist_PowerSpectrum)
parID_ps[parID] = ps
# print(final_concentration0)
if save==True:
if lhs==True:
filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)
pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )
# if root=='/rds/general/user/mo2016':
# sendmail('general_filename')
0%| | 0/20 [00:00<?, ?it/s]
/rds/general/user/mo2016
100%|██████████| 20/20 [00:03<00:00, 6.41it/s]
#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
modelling_home = '/Volumes/mo2016/home/Documents/modelling'
modelling_local = root + '/Documents/modelling'
import matplotlib as mpl
mpl.use('tkagg')
if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
modelling_ephemeral = root + '/ephemeral/Documents/modelling'
modelling_home = root + '/home/Documents/modelling'
modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)
from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################
%matplotlib inline
plot=True
save=True
if plot==True:
save==False
# %matplotlib inline
def fourier2d(final_concentration):
# print('Normalised PowerSpectrum of original')
pixels=final_concentration[5]
pixels[pixels == 0] = 0.0000000001
pixels/=np.sum(pixels) #line added
#fourier transform
f = np.fft.ifftshift(pixels)
f = np.fft.fft2(f)
f = np.fft.fftshift(f)
if plot==True:
plt.subplot(142)
plt.imshow(np.real(f))
plt.colorbar()
plt.subplot(143)
plt.imshow(np.imag(f))
plt.colorbar()
# max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
PowerSpectrum = np.real(f)**2 + np.imag(f)**2
# PowerSpectrum[max_coordinates]=0
if plot==True:
plt.subplot(144)
plt.imshow(PowerSpectrum)
plt.colorbar()
if np.sum(PowerSpectrum)!=0:
PowerSpectrum/=np.sum(PowerSpectrum)
# print('sum',np.sum(PowerSpectrum))
# print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))
# plt.subplot(132)
# plt.imshow(PowerSpectrum)#,norm=LogNorm())
# plt.colorbar()
# binwidth = 0.00001
PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
# print(np.sum(PowerSpectrum_vector))
# plt.subplot(133)
# hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
# plt.yscale('log')
# plt.show()
return PowerSpectrum_vector
def entropy(histogram):
x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
return x
lhs=True
if lhs==True:
folder = 'fullcircuit_5716gaussian/1M_turingI'
variant=0
else:
var=float(sys.argv[1])
# var=0.23
folder = 'fullcircuit_5716gaussian/var%s'%var
variant='5716gaussian'
circuit_n=2
shape='square'
mechanism = 'fullcircuit'
L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )
parID_ps = {}
# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
# print('parID',parID)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
else:
filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
if not np.isnan(final_concentration0).any():
final_concentration = np.round(final_concentration0,4)
if plot==True:
plt.figure(figsize=[24,4])
plt.subplot(141)
plt.imshow(final_concentration[5])
plt.colorbar()
hist_PowerSpectrum = fourier2d(final_concentration)
ps = entropy(hist_PowerSpectrum)
parID_ps[parID] = ps
# print(final_concentration0)
if save==True:
if lhs==True:
filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)
pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )
# if root=='/rds/general/user/mo2016':
# sendmail('general_filename')
5%|▌ | 1/20 [00:00<00:02, 7.40it/s]
/rds/general/user/mo2016
100%|██████████| 20/20 [00:02<00:00, 8.44it/s]
#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
modelling_home = '/Volumes/mo2016/home/Documents/modelling'
modelling_local = root + '/Documents/modelling'
import matplotlib as mpl
mpl.use('tkagg')
if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
modelling_ephemeral = root + '/ephemeral/Documents/modelling'
modelling_home = root + '/home/Documents/modelling'
modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)
from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################
%matplotlib inline
plot=True
save=True
if plot==True:
save==False
# %matplotlib inline
def fourier2d(final_concentration):
# print('Normalised PowerSpectrum of original')
pixels=final_concentration[5]
pixels[pixels == 0] = 0.0000000001
pixels/=np.sum(pixels) #line added
#fourier transform
f = np.fft.ifftshift(pixels)
f = np.fft.fft2(f)
f = np.fft.fftshift(f)
if plot==True:
plt.subplot(142)
plt.imshow(np.real(f))
plt.colorbar()
plt.subplot(143)
plt.imshow(np.imag(f))
plt.colorbar()
# max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
PowerSpectrum = np.real(f)**2 + np.imag(f)**2
# PowerSpectrum[max_coordinates]=0
# if plot==True:
# plt.subplot(144)
# plt.imshow(PowerSpectrum)
# plt.colorbar()
if np.sum(PowerSpectrum)!=0:
PowerSpectrum/=np.sum(PowerSpectrum)
if plot==True:
plt.subplot(144)
plt.imshow(PowerSpectrum)
plt.colorbar()
# print('sum',np.sum(PowerSpectrum))
# print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))
# plt.subplot(132)
# plt.imshow(PowerSpectrum)#,norm=LogNorm())
# plt.colorbar()
# binwidth = 0.00001
PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
# print(np.sum(PowerSpectrum_vector))
# plt.subplot(133)
# hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
# plt.yscale('log')
# plt.show()
return PowerSpectrum_vector
def entropy(histogram):
x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
return x
lhs=True
if lhs==True:
folder = 'fullcircuit_5716gaussian/1M_turingI'
variant=0
else:
var=float(sys.argv[1])
# var=0.23
folder = 'fullcircuit_5716gaussian/var%s'%var
variant='5716gaussian'
circuit_n=2
shape='square'
mechanism = 'fullcircuit'
L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )
parID_ps = {}
# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
# print('parID',parID)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
else:
filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
if not np.isnan(final_concentration0).any():
final_concentration = np.round(final_concentration0,4)
if plot==True:
plt.figure(figsize=[24,4])
plt.subplot(141)
plt.imshow(final_concentration[5])
plt.colorbar()
hist_PowerSpectrum = fourier2d(final_concentration)
ps = entropy(hist_PowerSpectrum)
parID_ps[parID] = ps
# print(final_concentration0)
if save==True:
if lhs==True:
filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)
pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )
# if root=='/rds/general/user/mo2016':
# sendmail('general_filename')
5%|▌ | 1/20 [00:00<00:01, 9.92it/s]
/rds/general/user/mo2016
100%|██████████| 20/20 [00:02<00:00, 7.57it/s]
#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
modelling_home = '/Volumes/mo2016/home/Documents/modelling'
modelling_local = root + '/Documents/modelling'
import matplotlib as mpl
mpl.use('tkagg')
if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
modelling_ephemeral = root + '/ephemeral/Documents/modelling'
modelling_home = root + '/home/Documents/modelling'
modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)
from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################
%matplotlib inline
plot=True
save=True
if plot==True:
save==False
# %matplotlib inline
def fourier2d(final_concentration):
# print('Normalised PowerSpectrum of original')
pixels=final_concentration[5]
pixels[pixels == 0] = 0.0000000001
pixels/=np.sum(pixels) #line added
#fourier transform
f = np.fft.ifftshift(pixels)
f = np.fft.fft2(f)
f = np.fft.fftshift(f)
if plot==True:
plt.subplot(142)
plt.imshow(np.real(f))
plt.colorbar()
plt.subplot(143)
plt.imshow(np.imag(f))
plt.colorbar()
max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
PowerSpectrum = np.real(f)**2 + np.imag(f)**2
PowerSpectrum[max_coordinates]=0
# if plot==True:
# plt.subplot(144)
# plt.imshow(PowerSpectrum)
# plt.colorbar()
if np.sum(PowerSpectrum)!=0:
PowerSpectrum/=np.sum(PowerSpectrum)
if plot==True:
plt.subplot(144)
plt.imshow(PowerSpectrum)
plt.colorbar()
# print('sum',np.sum(PowerSpectrum))
# print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))
# plt.subplot(132)
# plt.imshow(PowerSpectrum)#,norm=LogNorm())
# plt.colorbar()
# binwidth = 0.00001
PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
# print(np.sum(PowerSpectrum_vector))
# plt.subplot(133)
# hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
# plt.yscale('log')
# plt.show()
return PowerSpectrum_vector
def entropy(histogram):
x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
return x
lhs=True
if lhs==True:
folder = 'fullcircuit_5716gaussian/1M_turingI'
variant=0
else:
var=float(sys.argv[1])
# var=0.23
folder = 'fullcircuit_5716gaussian/var%s'%var
variant='5716gaussian'
circuit_n=2
shape='square'
mechanism = 'fullcircuit'
L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )
parID_ps = {}
# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
# print('parID',parID)
if lhs==True:
filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
else:
filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
if not np.isnan(final_concentration0).any():
final_concentration = np.round(final_concentration0,4)
if plot==True:
plt.figure(figsize=[24,4])
plt.subplot(141)
plt.imshow(final_concentration[5])
plt.colorbar()
hist_PowerSpectrum = fourier2d(final_concentration)
ps = entropy(hist_PowerSpectrum)
parID_ps[parID] = ps
# print(final_concentration0)
if save==True:
if lhs==True:
filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)
pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )
# if root=='/rds/general/user/mo2016':
# sendmail('general_filename')
5%|▌ | 1/20 [00:00<00:01, 9.63it/s]
/rds/general/user/mo2016
100%|██████████| 20/20 [00:02<00:00, 7.95it/s]